########################################################################
# PLS 206 — Week 3: General(ized) linear models
# Grey Monroe, Fall 2026
#
# Every block is tagged:
#   [CORE]   the analysis: code you would run on your own data
#   [HANDY]  useful, not essential: wrangling, plotting, small helpers
#   [DEMO]   for the lecture only: simulated data, sums by hand, and the
#            code behind slide figures
########################################################################

# [CORE]
# install.packages(c("ggplot2", "emmeans", "car"))
library(ggplot2)

########################################################################
# PART A. Groups as predictors
########################################################################

# ---- 1) The data: a wheat field trial --------------------------------

# 90 plots: three varieties (Amber, Bronze, Copper) x two water
# treatments (control, drought), 15 plots of each, and a nitrogen rate
# (kg N/ha) that varies from plot to plot. Three responses: yield
# (t/ha), rust (stripe rust seen on the plot: 1 = yes, 0 = no), and
# aphids (aphids counted on 10 tillers). Simulated.

# [CORE]
trial <- read.csv("https://greymonroe.github.io/PLS_206/data/wheat_trial.csv")
head(trial)
table(trial$variety, trial$water)   # check the design: 15 plots per combination

# [HANDY]
# groups: violin + box + points
ggplot(trial, aes(water, yield, fill = water, colour = water)) +
  geom_violin(alpha = 0.25, colour = NA) +
  geom_boxplot(width = 0.18, fill = "white", colour = "grey30", outlier.shape = NA) +
  geom_jitter(width = 0.08, height = 0) +
  labs(x = NULL, y = "Yield (t/ha)") +
  theme_classic() + theme(legend.position = "none")

# ---- 2) Two groups: a dummy variable ---------------------------------

# A group becomes a 0/1 column: 0 for control, 1 for drought. The
# intercept is the control mean; the coefficient is drought minus
# control.

# [CORE]
m_water <- lm(yield ~ water, data = trial)
round(coef(summary(m_water)), 3)
tapply(trial$yield, trial$water, mean)   # the two means

# [DEMO: the drought mean from the coefficients]
coef(m_water)[1] + coef(m_water)[2]   # intercept + difference = drought mean

# ---- 3) It's a t-test ------------------------------------------------

# [CORE]
t.test(yield ~ water, data = trial, var.equal = TRUE)   # same t and p as waterdrought

# ---- 4) Under the hood: the design matrix ----------------------------

# [CORE]
head(trial$water)
head(model.matrix(m_water))   # the 0/1 column R made for you

# ---- 5) Which group is the intercept? --------------------------------

# Text columns go in alphabetical order; the first level is the
# reference.

# [CORE]
trial2 <- trial
trial2$water <- factor(trial2$water,
                       levels = c("drought", "control"))
coef(lm(yield ~ water, data = trial2))   # same model, flipped sign

########################################################################
# PART B. More than two groups (ANOVA)
########################################################################

# ---- 6) Three varieties, two dummies ---------------------------------

# The 45 control plots only, for now.

# [CORE]
control <- subset(trial, water == "control")
tapply(control$yield, control$variety, mean)
m_var <- lm(yield ~ variety, data = control)
round(coef(summary(m_var)), 3)   # each row: that variety minus Amber

# [DEMO: the design matrix with three groups]
mm <- model.matrix(m_var)
head(cbind(control["variety"], mm), 7)   # Amber plots: zeros in both dummies

# ---- 7) Do the varieties differ at all? One F-test -------------------

# [CORE]
anova(m_var)                     # the one-way ANOVA
summary(aov(yield ~ variety, data = control))   # the same table, Week 1 style

# ---- 8) Which pairs? Tukey, and adjusting for many tests -------------

# [CORE]
TukeyHSD(aov(m_var))

# [HANDY]
# raw vs adjusted p-values for every pair
pairwise.t.test(control$yield, control$variety, p.adjust.method = "none")
pairwise.t.test(control$yield, control$variety, p.adjust.method = "holm")
# p.adjust() works on any vector of p-values
p <- c(0.004, 0.02, 0.03, 0.2)    # any set of p-values
p.adjust(p, method = "holm")
p.adjust(p, method = "BH")         # false discovery rate

# [DEMO: the chance of at least one false positive in k tests]
k <- 1:20
plot(k, 1 - 0.95^k, type = "b", col = "red", xlab = "Number of tests",
     ylab = "Chance of at least one false positive")
abline(h = 0.05, lty = 2)

# ---- 9) How big? Eta-squared -----------------------------------------

# [CORE]
ss <- anova(m_var)$"Sum Sq"
ss[1] / sum(ss)              # eta-squared: share of variation explained by variety
summary(m_var)$r.squared     # the same, with one factor

########################################################################
# PART C. Two factors
########################################################################

# ---- 10) Additive: one drought effect for every variety --------------

# [HANDY]
ggplot(trial, aes(variety, yield, colour = water)) +
  geom_point(position = position_jitterdodge(jitter.width = 0.15, dodge.width = 0.6)) +
  labs(x = NULL, y = "Yield (t/ha)") +
  theme_classic()

# [CORE]
m_add <- lm(yield ~ variety + water, data = trial)
round(coef(summary(m_add)), 3)

# ---- 11) Interaction: each variety gets its own drought effect -------

# [CORE]
m_int <- lm(yield ~ variety * water, data = trial)
round(coef(summary(m_int)), 3)
anova(m_int)   # test the interaction row first

# [DEMO: a cell mean from the coefficients]
b <- coef(m_int)
b["(Intercept)"] + b["varietyCopper"] + b["waterdrought"] + b["varietyCopper:waterdrought"]   # Copper, drought
mean(trial$yield[trial$variety == "Copper" & trial$water == "drought"])

# [HANDY]
# interaction plot: predicted means, one line per variety
g <- expand.grid(variety = c("Amber", "Bronze", "Copper"), water = c("control", "drought"))
g$additive    <- predict(m_add, g)
g$interaction <- predict(m_int, g)
ggplot(g, aes(water, interaction, colour = variety, group = variety)) +
  geom_line(linewidth = 1.2) + geom_point(size = 3) +
  labs(x = NULL, y = "Predicted yield (t/ha)") +
  theme_classic()

# ---- 11b) summary() vs anova(): what is an "effect"? -----------------

# Every row of summary() is one comparison between two cells. Every row
# of anova() is a whole term. With an interaction in the model, the
# "main effect" coefficients are effects at the reference level of the
# other factor: varietyCopper is Copper vs Amber under control only,
# and waterdrought is the drought effect in Amber only.

# [CORE]
round(coef(summary(m_int)), 3)[, c(1, 4)]   # single comparisons
anova(m_int)                                # whole terms

# Change the reference level: the coefficient changes (even its sign);
# the ANOVA table does not.

# [CORE]
trial2$water <- factor(trial$water, levels = c("drought", "control"))
m_int2 <- lm(yield ~ variety * water, data = trial2)
round(coef(summary(m_int))["varietyCopper", ], 3)    # Copper - Amber, under control
round(coef(summary(m_int2))["varietyCopper", ], 3)   # Copper - Amber, under drought
anova(m_int)["variety", "Pr(>F)"]
anova(m_int2)["variety", "Pr(>F)"]                   # identical

# The anova() variety row is variety averaged over water. With a real
# interaction, that average can hide the story (emmeans warns you):

# [CORE]
emmeans::emmeans(m_int, ~ variety)   # averaged over water

# ---- 12) The effect within each variety: emmeans ---------------------

# [CORE]
library(emmeans)
emmeans(m_int, ~ water | variety)          # model-predicted means with CIs
pairs(emmeans(m_int, ~ water | variety))   # the drought effect in each variety

# ---- 13) A second factor removes noise -------------------------------

# [CORE]
anova(lm(yield ~ variety, data = trial))   # drought left in the residuals
anova(m_add)                               # drought accounted for

# ---- 14) Unequal group sizes: order matters --------------------------

# [CORE]
unbal <- trial[-(1:20), ]        # lose 20 plots
table(unbal$variety, unbal$water)
anova(lm(yield ~ variety + water, data = unbal))   # sequential (Type I)
anova(lm(yield ~ water + variety, data = unbal))   # the other order
library(car)
Anova(lm(yield ~ variety + water, data = unbal), type = 2)   # each adjusted for the other

# ---- 14a) Overlapping factors: who gets the shared slice? ------------

# A lopsided subset of the trial, like on-farm data: drop 12 of the 15
# Bronze drought plots and 12 of the 15 Copper control plots, so Bronze
# is mostly irrigated and Copper mostly droughted. Now variety and
# water overlap, and part of the variation could be credited to either.
# Type I gives that shared slice to whichever term is first; Type II
# gives it to neither.

# [DEMO: a lopsided subset of the trial]
drop <- c(which(trial$variety == "Bronze" & trial$water == "drought")[-(1:3)],
          which(trial$variety == "Copper" & trial$water == "control")[-(1:3)])
lopsided <- trial[-drop, ]

# [CORE]
table(lopsided$variety, lopsided$water)
anova(lm(yield ~ variety + water, data = lopsided))   # variety first: gets the shared slice
anova(lm(yield ~ water + variety, data = lopsided))   # water first: variety keeps only its own
Anova(lm(yield ~ variety + water, data = lopsided), type = 2)   # each on its own slice

# ---- 14b) Type III, and the contrasts trap ---------------------------

# Type III adjusts every term for every other term, the interaction
# included. With an interaction, variety has a different effect at each
# water level, so Type III has to test it somewhere. With R's default
# contrasts, its "variety" row tests variety at the reference water
# level (like the summary() row), so the answer changes with the
# reference. With sum-to-zero contrasts it tests variety averaged over
# water. This happens in balanced data too: it comes from the
# interaction, not from unequal group sizes.

# [CORE]
library(car)
Anova(m_int, type = 3)["variety", 3:4]    # control as reference: F = 58
Anova(m_int2, type = 3)["variety", 3:4]   # drought as reference: nothing
options(contrasts = c("contr.sum", "contr.poly"))
m_sum <- lm(yield ~ variety * water, data = trial)
Anova(m_sum, type = 3)                    # sum-to-zero: variety averaged over water
options(contrasts = c("contr.treatment", "contr.poly"))   # back to the default
emmeans::joint_tests(m_int)               # the same averaged tests, whatever the coding

# [DEMO: what the variety columns are under each coding]
cells <- data.frame(variety = rep(c("Amber", "Bronze", "Copper"), 2), water = rep(c("control", "drought"), each = 3))
model.matrix(~ variety * water, cells)        # treatment: interaction columns are 0 for every control row
options(contrasts = c("contr.sum", "contr.poly"))
model.matrix(~ variety * water, cells)        # sum-to-zero: +1/-1 in both water levels
options(contrasts = c("contr.treatment", "contr.poly"))

# Which type? Balanced design: Types I, II, and III (with contr.sum)
# agree. Unbalanced: Type I only if the order means something; Type II
# when there is no interaction; Type III with contr.sum if you must
# test main effects alongside an interaction. If the interaction is
# real, report the effects within each level (emmeans) instead.

########################################################################
# PART D. Groups and a slope together (ANCOVA)
########################################################################

# ---- 15) Parallel lines ----------------------------------------------

# [CORE]
m_anc <- lm(yield ~ nitrogen + water, data = trial)
round(coef(summary(m_anc)), 4)   # intercept = control at 0 kg N

# ---- 16) Different slopes --------------------------------------------

# [CORE]
m_anc2 <- lm(yield ~ nitrogen * water, data = trial)
round(coef(summary(m_anc2)), 4)
# nitrogen slope under drought: nitrogen + nitrogen:waterdrought
coef(m_anc2)["nitrogen"] + coef(m_anc2)["nitrogen:waterdrought"]

# With the interaction in, waterdrought is the drought effect at 0 kg
# N. Centre nitrogen and it becomes the drought effect at the average
# rate:

# [CORE]
trial$nitrogen_c <- trial$nitrogen - mean(trial$nitrogen)
m_anc3 <- lm(yield ~ nitrogen_c * water, data = trial)
round(coef(summary(m_anc3)), 4)

# [HANDY]
# one fitted line per group
ggplot(trial, aes(nitrogen, yield, colour = water)) +
  geom_point() +
  geom_smooth(method = "lm", se = FALSE) +   # separate slopes, like nitrogen * water
  labs(x = "Nitrogen (kg N/ha)", y = "Yield (t/ha)") +
  theme_classic()

# ---- 16b) Two numeric predictors: the main-effect t-test is the slope at 0 ---

# With nitrogen * water, the nitrogen coefficient is the nitrogen slope
# where water = 0. Here no plot is that dry, and the coefficient looks
# like nothing while anova() calls nitrogen one of the strongest terms.

# [DEMO: simulate an irrigation x nitrogen trial]
set.seed(47)
irrig <- data.frame(nitrogen = round(runif(60, 0, 200)),
                    water = round(runif(60, 15, 45)))   # mm/week
irrig$yield <- round(2 + 0.04 * irrig$water +
  0.0004 * irrig$nitrogen * irrig$water + rnorm(60, 0, 0.6), 2)

# [CORE]
m_nw <- lm(yield ~ nitrogen * water, data = irrig)
round(coef(summary(m_nw)), 4)   # nitrogen: the slope at water = 0
anova(m_nw)                     # nitrogen: a whole term, entered first

# centre both, refit: same fit, same interaction; main effects now at the means
irrig$nitrogen_c <- irrig$nitrogen - mean(irrig$nitrogen)
irrig$water_c    <- irrig$water - mean(irrig$water)
m_nwc <- lm(yield ~ nitrogen_c * water_c, data = irrig)
round(coef(summary(m_nwc)), 4)
c(summary(m_nw)$r.squared, summary(m_nwc)$r.squared)   # identical

# [DEMO: the nitrogen slope at every water level, and how the tests line up]
b <- coef(m_nw); V <- vcov(m_nw)
w <- seq(0, 50, by = 0.5)
slope <- b["nitrogen"] + b["nitrogen:water"] * w
se <- sqrt(V["nitrogen", "nitrogen"] + w^2 * V["nitrogen:water", "nitrogen:water"] + 2 * w * V["nitrogen", "nitrogen:water"])
plot(w, slope, type = "l", lwd = 2, col = "blue", ylim = range(slope - 2 * se, slope + 2 * se),
     xlab = "Water (mm/week)", ylab = "Nitrogen slope (t/ha per kg N)")
lines(w, slope - 2 * se, lty = 3, col = "blue"); lines(w, slope + 2 * se, lty = 3, col = "blue")
abline(h = 0, lty = 2); abline(v = range(irrig$water), col = "grey")
# centred t^2 = Type III F; Type I and II test nitrogen without the interaction
coef(summary(m_nwc))["nitrogen_c", "t value"]^2
library(car)
Anova(m_nwc, type = 3)["nitrogen_c", "F value"]
Anova(m_nw, type = 2)["nitrogen", "F value"]
Anova(m_nw, type = 3)["nitrogen", "F value"]   # uncentred Type III: the same at-zero test

########################################################################
# PART E. Yes/no responses: logistic regression
########################################################################

# ---- 17) A straight line fails ---------------------------------------

# [CORE]
m_lin <- lm(rust ~ nitrogen + variety + water, data = trial)
range(fitted(m_lin))                          # fitted "probabilities" for the real plots
sum(fitted(m_lin) < 0 | fitted(m_lin) > 1)   # how many are impossible

# ---- 18) Probability, odds, log-odds ---------------------------------

# [CORE]
p <- c(0.2, 0.5, 0.8, 0.99)
odds <- p / (1 - p)
cbind(p, odds, log_odds = log(odds))
qlogis(0.8)    # probability -> log-odds
plogis(1.39)   # log-odds -> probability

# [DEMO: the logit curve and the S-curve]
par(mfrow = c(1, 2))
curve(qlogis(x), 0.005, 0.995, xlab = "Probability", ylab = "Log-odds", col = "blue", lwd = 2)
curve(plogis(x), -6, 6, xlab = "Log-odds", ylab = "Probability", col = "red", lwd = 2)
par(mfrow = c(1, 1))

# ---- 19) A group as the predictor ------------------------------------

# [CORE]
table(trial$rust, trial$water)
m_rw <- glm(rust ~ water, family = binomial, data = trial)
round(coef(summary(m_rw)), 3)   # log-odds scale
plogis(coef(m_rw)[1])           # P(rust), control
plogis(sum(coef(m_rw)))         # P(rust), drought
prop.table(table(trial$rust, trial$water), 2)
exp(coef(m_rw))                 # odds ratio for drought
exp(confint(m_rw))

# [HANDY]
# the older tools for a 2 x 2 table ask the same question
fisher.test(table(trial$water, trial$rust))
chisq.test(table(trial$water, trial$rust))

# ---- 20) A continuous predictor --------------------------------------

# [CORE]
m_rn <- glm(rust ~ nitrogen, family = binomial, data = trial)
round(coef(summary(m_rn)), 4)
exp(10 * coef(m_rn)["nitrogen"])   # odds ratio per 10 kg N
-coef(m_rn)[1] / coef(m_rn)[2]     # nitrogen where P(rust) = 0.5

# [HANDY]
# plot the fitted curve: predict probabilities over a grid
grid <- data.frame(nitrogen = 0:200)
grid$p <- predict(m_rn, newdata = grid, type = "response")
plot(trial$nitrogen, jitter(trial$rust, 0.1), xlab = "Nitrogen (kg N/ha)", ylab = "Rust (0/1)")
lines(grid$nitrogen, grid$p, col = "red", lwd = 2)

# ---- 21) predict(): log-odds or probability --------------------------

# [CORE]
new <- data.frame(nitrogen = c(50, 100, 150))
predict(m_rn, newdata = new)                     # log-odds (the default)
predict(m_rn, newdata = new, type = "response")  # probability
head(fitted(m_rn))                               # fitted() is always probability

# ---- 22) How often is it right? --------------------------------------

# [CORE]
pred <- as.numeric(fitted(m_rn) > 0.5)
table(predicted = pred, observed = trial$rust)   # confusion matrix
mean(pred == trial$rust)    # accuracy
mean(trial$rust == 0)       # always guess "no rust": the baseline to beat

# ---- 23) Several predictors ------------------------------------------

# [CORE]
m_rall <- glm(rust ~ nitrogen + variety + water, family = binomial, data = trial)
round(coef(summary(m_rall)), 3)
round(exp(coef(m_rall)), 2)   # odds ratios, each holding the others fixed

# ---- 24) Deviance ----------------------------------------------------

# [CORE]
m_rall$null.deviance                  # intercept only
deviance(m_rall)                      # with the predictors
1 - deviance(m_rall) / m_rall$null.deviance   # McFadden's pseudo-R2
anova(m_rn, m_rall, test = "Chisq")   # nested GLMs: likelihood-ratio test

# ---- 25) Counts out of a total: cbind(successes, failures) -----------

# [DEMO: simulate germination in 30 trays of 20 seeds]
set.seed(3)
trays <- data.frame(temp = runif(30, 8, 30))
trays$germ <- rbinom(30, 20, plogis(-4 + 0.25 * trays$temp))

# [CORE]
m_g <- glm(cbind(germ, 20 - germ) ~ temp, family = binomial, data = trays)
round(coef(summary(m_g)), 3)
exp(coef(m_g)["temp"])   # odds of germinating, per degree

# ---- 26) More than two categories: multinomial -----------------------

# [HANDY]
library(nnet)
m_sp <- multinom(Species ~ Petal.Length + Petal.Width, data = iris, trace = FALSE)
table(predicted = predict(m_sp), observed = iris$Species)

########################################################################
# PART F. Counts: Poisson regression
########################################################################

# ---- 27) Last week's seed counts -------------------------------------

# [DEMO: simulate last week's seed counts]
set.seed(4)
size  <- runif(80, 1, 10)
seeds <- rpois(80, lambda = exp(-1 + 0.35 * size))
pods  <- data.frame(size, seeds)

# [CORE]
m_pod <- glm(seeds ~ size, family = poisson, data = pods)
round(coef(summary(m_pod)), 3)   # log scale
exp(coef(m_pod))                 # rate ratio: x 1.46 seeds per unit of size

# [HANDY]
# the Poisson curve vs last week's straight line
grid <- data.frame(size = seq(0, 10, by = 0.1))
plot(pods$size, pods$seeds, xlab = "Plant size", ylab = "Seeds per plant")
abline(lm(seeds ~ size, data = pods), lty = 2)
lines(grid$size, predict(m_pod, newdata = grid, type = "response"), col = "red", lwd = 2)

# ---- 28) Aphids, and overdispersion ----------------------------------

# [CORE]
m_aph <- glm(aphids ~ nitrogen, family = poisson, data = trial)
round(coef(summary(m_aph)), 4)
deviance(m_pod) / df.residual(m_pod)   # about 1: fine
deviance(m_aph) / df.residual(m_aph)   # about 6: overdispersed

# ---- 29) Fix: negative binomial (or quasipoisson) --------------------

# [CORE]
library(MASS)
m_nb <- glm.nb(aphids ~ nitrogen, data = trial)
round(coef(summary(m_nb)), 4)   # same slope, honest SE

# [HANDY]
m_qp <- glm(aphids ~ nitrogen, family = quasipoisson, data = trial)
round(coef(summary(m_qp)), 4)   # same slope as Poisson, SE scaled by the dispersion

########################################################################
# PART G. Reporting
########################################################################

# ---- 30) Odds ratios and rate ratios with CIs, on a sensible step ----

# [CORE]
exp(10 * coef(m_rn)["nitrogen"])
exp(10 * confint(m_rn)["nitrogen", ])
# "The odds of rust increased with nitrogen (odds ratio 1.16 per 10 kg N,
#  95% CI 1.07 to 1.28; z = 3.5, p < 0.001; logistic regression, n = 90)."
exp(10 * coef(m_nb)["nitrogen"])
exp(10 * confint(m_nb)["nitrogen", ])
# "Aphid counts increased with nitrogen (rate ratio 1.11 per 10 kg N,
#  95% CI 1.07 to 1.14; negative binomial GLM, n = 90)."

########################################################################
# PART H. Bonus: mixed effect models
########################################################################

# ---- 31) Plots grouped on farms --------------------------------------

# 112 plots on 12 farms (2 to 25 plots each). Nitrogen raises yield
# everywhere, each farm has its own baseline, and cover crops (chosen
# per farm) do nothing.

# [DEMO: simulate 112 plots on 12 farms]
set.seed(181)
nj <- c(2, 2, 3, 4, 5, 6, 8, 10, 12, 15, 20, 25)   # plots per farm
ids <- sprintf("farm%02d", 1:12)
cover_j <- setNames(rep(c("no", "yes"), 6)[sample(12)], ids)
u <- setNames(rnorm(12, 0, 0.8), ids)              # each farm: its own baseline
farm <- rep(ids, nj)
nitrogen <- round(runif(112, 0, 200))
yield <- round(3.5 + 0.012 * nitrogen + u[farm] + rnorm(112, 0, 1), 2)
farms <- data.frame(farm, cover = unname(cover_j[farm]), nitrogen, yield)

# [CORE]
head(farms, 4)
table(farms$farm)

# ---- 32) Farm as fixed vs random -------------------------------------

# [CORE]
# install.packages("lme4")
library(lme4)
m_fixed <- lm(yield ~ nitrogen + farm, data = farms)          # 12 free intercepts
m_mixed <- lmer(yield ~ nitrogen + (1 | farm), data = farms)  # farms drawn from N(0, sd_farm^2)
VarCorr(m_mixed)                       # two variance components
round(coef(summary(m_mixed)), 4)       # fixed effects (no p-values: see lmerTest)
round(coef(summary(m_fixed))["nitrogen", ], 4)

# ---- 33) Shrinkage: small farms borrow strength ----------------------

# [DEMO: farm effects, fixed vs random]
ids <- sort(unique(farms$farm))
fixed  <- c(0, coef(m_fixed)[paste0("farm", ids[-1])]) + coef(m_fixed)[1] - fixef(m_mixed)[1]
random <- ranef(m_mixed)$farm[ids, 1]
n <- as.numeric(table(farms$farm)[ids])
plot(1:12, fixed, ylim = range(c(fixed, random)), xaxt = "n", xlab = "", ylab = "Farm effect (t/ha)")
axis(1, at = 1:12, labels = paste0("n=", n), cex.axis = 0.7)
points(1:12, random, pch = 16, col = "purple")
arrows(1:12, fixed, 1:12, random, length = 0.05, col = "purple")
abline(h = 0, lty = 2)
# the weight each farm keeps on its own estimate
sd_farm <- attr(VarCorr(m_mixed)$farm, "stddev"); s <- sigma(m_mixed)
round(sd_farm^2 / (sd_farm^2 + s^2 / n), 2)

# ---- 34) Why it matters: a treatment chosen per farm -----------------

# [CORE]
round(coef(summary(lm(yield ~ nitrogen + cover, data = farms)))["coveryes", ], 4)   # plots as independent: wrong
m_mixed2 <- lmer(yield ~ nitrogen + cover + (1 | farm), data = farms)
round(coef(summary(m_mixed2))["coveryes", ], 4)                                     # farm in the model: honest

# [HANDY]
# checks: residuals as for lm(), and the farm effects should look normal
plot(m_mixed)
qqnorm(ranef(m_mixed)$farm[, 1]); qqline(ranef(m_mixed)$farm[, 1])
# the same for yes/no or counts: glmer(rust ~ nitrogen + (1 | block), family = binomial, data = ...)